#include <iostream>
#include <cstdlib>
#include <ctime>
#include <cmath>

using namespace std;


#define cutoffw 1.0
double stepsize = 0.5;
double w = 0.01;
long acceptMove = 0;
double *data = NULL;


double getProbability(double x){
    
    if (x < 0){
        return 0;
    }
    else{
        return x * exp(-1.0 *  x / cutoffw);
    }

}

double ran1(){
    return 1.0 * rand() / RAND_MAX;
}

double ran(){
    return (ran1() - 0.5);
}

double runToEquilibrium(long relax){

    for (long i = 0; i < relax; i++){
        double trial = w + ran() * stepsize;
        double ratio = getProbability(trial) / getProbability(w);
        if (ratio > 1 || ratio > ran1()){
            w = trial;
            acceptMove++;
        }
    }

    double temp = 1.0 * acceptMove / relax;
    cout << "relax = " << relax << " stepsize = " << stepsize << "\tacceptRatio = " << temp << endl;
    acceptMove = 0;
    return temp;
}


double generateDistribution(long num = 20){
    data = new double[20];
    for (long i = 0; i < num; i++){
        data[i] = 0;
    }

    for(long i = 0; i < num; i++){
        double trial = w + ran() * stepsize; 
        double ratio = getProbability(trial) / getProbability(w);
        if (ratio > 1 || ratio > ran1()){
            w = trial;
            acceptMove++;
        }
        data[i] = w;
    }
}


void print(long num){
    for (long i  = 0; i < num; i++){
        cout << data[i] <<  " ";
    }
    cout << endl;
}




int main(){
    
    srand(time(NULL));
    double temp = runToEquilibrium(20000);
    for (;;){
        if (temp > 0.48){
            stepsize = stepsize * 1.05;
            temp = runToEquilibrium(20000);
        }
        else if (temp < 0.4){
            stepsize = stepsize * 0.95;
            temp = runToEquilibrium(20000);
        }
        else{
            break;
        }
    }
    generateDistribution(20);
    print(20);
    delete []data;
    return 0;
}







